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1.0  INTRODUCTION 


A  study  titled  "Imagery  Exploitation  Using  Wavelets"  was  conducted  at 
Westinghouse  Electric  Corporation  to  further  the  exploitation  of  digital 
reconnaissance  imagery.  The  study  results  provide  a  set  of  wavelet  based  tools  to  be 
included  in  the  Image  Processing  Toolkit.  The  tools  were  designed  to  aid  in 
visualization  of  the  interaction  of  the  wavelet  transform  and  imagery,  to  provide 
different  views  of  the  imagery  (such  as  pseudocoloration)  and  to  enhance  certain 
image  features  based  on  the  localization  property  of  the  wavelet  transform.  This 
Final  Report  provides  some  lessons  learned  during  the  development  of  the  tools,  a 
general  description  of  the  tools  and  recommendations  for  further  work  in  this 
area. 


2.0  OBJECTIVE 

The  objective  of  this  study  was  to  develop  a  set  of  image  exploitation  tools  using 
wavelets.  Earlier  experimental  work  suggested  that  tailoring  the  mother  wavelet  to 
the  task  at  hand  could  enhance  the  ability  to  extract  certain  image  features.  For 
example,  this  was  shown  in  the  extraction  of  road  scarring  in  Synthetic  Appature 
RADAR  (SAR)  imagery.  In  view  of  the  above,  two  classes  of  wavelets  were 
implemented  in  the  study:  general  frame  wavelets  and  orthonormal  wavelets.  The 
wavelets,  which  are  covered  in  the  next  section,  are  applicable  to  most  of  the  tools. 
Six  tools  were  developed  and  are  listed  below: 

WAV_SCA  -  Wavelet  Scalogram 

WAV_CEP  -  Wavelet  Cepstrum  (Originally  named  WAV_CIC) 

WAV_MAC  -  Wavelet  Multiresolution  Automatic  Color  Rendering 
WAV_VOL  -  Wavelet  Volumepass  Filtering 
WAV_TEN  -  Wavelet  Transform  Encoding 
WAV_SOF  -  Wavelet  Scale  Orientation  Filter 

3.0  WAVELET  FUNCTION  AND  IMPLEMENTATION 

In  a  word,  wavelets  are  locally  supported  functions  labeled  with  scale  and  delay 
and  comprising  a  basis  for  function  space.  First  a  mother  wavelet  (Figure  3.0-1) 
is  selected  according  to  an  admissibility  condition;  then,  the  mother  is  labeled 
with  scale  and  delay  to  generate  the  wavelets  which  comprise  a  basis  (Figure 
3.0-2). 

Wavelets  are  attractive  for  the  exploitation  of  digital  reconnaissance  imagery  for 
many  reasons.  Wavelets  may  be  used  to  discern  the  local  frequency  content  of  a 
signal.  Data  may  be  decomposed  into  its  scale  and  spatial-orientation 
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Figure  3.0-1  Example  Mother  Wavelets  h(t) 


components  using  wavelets  [3].  The  mathematical  rigor  of  the  wavelets 
admissibility  and  frame  condition  precludes  the  inadvertent  loss  of  information  or 
the  introduction  of  artifacts.  The  frame  condition  [1]  equips  both  general-frame 
and  orthonormal  wavelets  with  fast  transforms  for  near-real  time  operation  and 
efficient  data  storage. 

By  far  the  most  important  attribute  of  wavelets,  however,  is  the  latitude  they  offer  to 
tailor  the  analysis  basis  to  suit  the  application  at  hand.  This  is  due  to  the  rather  lax 
admissibility  condition  on  the  mother  wavelet  [2]  which  requires  that  the  inner 
product  between  the  inverse  magnitude  of  frequency  and  the  magnitude  squared  of 
the  Fourier  transform  of  the  mother  be  finite: 

(  -i  A  o 

\  |  CO  |  I  h(co)  |  dco  <  OO  Eq.  3.0.1  Mother  Admissibility 

This  condition  is  satisfied  by  an  infinite  variety  of  functions.  Figure  3.0-1 
highlights  the  diversity  among  allowable  mother  wavelets.  This  freedom  to  tailor 
the  wavelets  to  suit  the  application  at  hand  makes  it  possible  to  invoke  wavelet  basis 
which,  during  transform,  force  application-important  features  to  congeal  in 
localized  regions  of  influence  in  the  wavelet  domain.  For  image  compression, 
coefficients  outside  the  region  of  influence  may  be  quantized  to  zero  while,  for 
image  exploitation,  coefficients  in  the  region  of  influence  may  be  filtered  using 
volumepass  or  cepstrum  filters.  However,  in  all  of  these  applications,  the  common 
practice  of  tailoring  the  mother  to  cause  qualitatively  important  data  features  to 
congeal  in  a  localized  region  in  the  wavelet  domain  is  exhibited. 

An  admissible  mother  h(t)  is  converted  to  a  basis  of  wavelets  by  substituting 
t — >(t-T)/a  in  h(t)  and  multiplying  by  lA/cT  for  normalization: 

1  t  - 1 

h  „  _(t)  =  ——  h( -  )  Eq.  3.0.2  Basis  of  Wavelets 

a’  Va  a 

Scale  a  being  in  the  denominator  in  the  argument  provides  the  contractions  and 
dilations  of  the  mother  while  %  being  subtracted  in  the  numerator  in  the  argument 
provides  time-shifted  versions  at  each  scale.  A  wavelet  transform  of  an  input 
function  x(t)  is  accomplished  by  taking  the  inner  product  between  the  input 
function  (or  sequence)  x(t)  and  all  of  the  labeled  wavelets  and  is  a  function  of  scale  a 
and  delay  t: 

Wavelet  Transform  (x(t))[a,r]  =  <  x(t),  ha  T(t)  >  Eq.  3.0.3  Wavelet  Transform 

While  the  mother  wavelet  must  satisfy  an  admissibility  condition,  the  labeled 
wavelets  themselves  must  meet  a  frame  condition  [1],  which  governs  the  spacing  in 
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the  aj-grid,  in  order  for  the  wavelet  transform  to  be  reversible  and 
information-preserving.  A  generalized  Parseval’s  Theorem,  it  requires  a 
proportionality  for  each  x(t)  between  the  energy  of  a  function  in  the  wavelets 
domain  and  its  norm  in  the  input  domain: 

A  ||x(t)||2  A  2  |<x(t),h  (t)  >|2  A  B  ||x(t)||2 

a,x  cc,t 

A  >  0,  B  <  OO  Eq.  3.0.4  Frame  Condition 

A  fast  wavelet  transform  (FWT)  may  be  realized  for  both  general-frame  and 
orthonormal  bases  of  wavelets  by  setting  the  a,t-grid  at  the  widest  possible 
spacing  (critical  frame).  This  process  is  illustrated  in  Figure  3.0-3  where  a 
square  wave  is  forward  -  and  then  inverse  -  transformed  for  three  a,i-spacings. 

To  set  a  wavelet  transform  for  maximum  speed  and  efficiency,  the  a,t-spacing 
(the  dots  on  the  left  of  Figure  3.0-3)  is  gradually  "tightened"  until  the 
reconstruction  (the  function  x(t)  on  the  right  of  Figure  3.0-3)  converges.  The 
exact  a,T-grid  spacing,  which  is  as  sparse  as  possible  while  maintaining 
reversibility,  is  called  critical  frame.  All  wavelet  transforms  (except  those  for 
scalograms)  shall  be  set  at  critical  frame. 


Figure  3.0-3  Setting  a  Wavelet  Transform  at  Critical  Frame. 


3.1  General  Frame  Wavelets 


The  following  General  Frame  Wavelets  meet  the  frame  and  admissibility 
conditions  and  were  implemented  in  this  study  : 


Haar  -  The  simplest  of  wavelets  assuming  the  values  of  plus  1  and  minus  1  over 

small  consecutive  intervals.  .  1 

1,  0  4  t  <  — 

-1.4-4  *  <  i 


Prototype: 


h  (t)  = 


{ 


0 ,  otherwise 


2SOPS  -  A  derivative  of  the  Haar  wavelet  using  a  second  order  low  pass  Filter  to 
define  the  wavelet  shape. 

t 

Prototype:  h  (t)  =  — - - -r 

(  t2  +  l)2 


Morlet  -  Sine  and  cosine  functions  which  are  windowed  by  a  Gaussian  to  limit 
their  extent.  The  transform  is  typically  implemented  with  complex  numbers. 


2nd  Dog  -  Second  derivative  of  a  Gaussian. 

_t_2 

Prototype:  h(t)  =  (l-t?)  e  2 


3rd  Dog  -  Third  derivative  of  a  Gaussian. 

_  t; 

Prototype:  h  (t)  =  ( t3  -  3t )  e  2 


Additional  developments  on  2SOPS  and  the  derivative  of  Gaussian  Wavelets  can 
be  found  in  document  "New  Classes  of  Frame  Wavelets  for  Applications"[4], 


With  a  change  in  scale,  the  wavelet  function  approaches  the  input  sample 
number  in  length.  As  a  consequence,  the  number  of  operations  involved  in 
computing  the  inner  products  becomes  very  large  and  results  in  long  processing 
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times  for  moderate  size  images.  Two  approaches  were  taken  in  order  to  speed  up 
the  process: 

1.  The  first  method  employed  the  critical  frame  process  discussed  previously. 
For  this  method  some  of  the  mother  wavelets  were  pre-computed  at  the 
various  non-integer  scale  and  delay  values  and  used  with  table  lookup 
procedures. 

2.  The  second  method  performed  the  correlations  using  FFT  methods  as 
described  below.  The  correlation  between  two  function  g  and  h  is  given  by: 

Corr(g,h)  =  ^g(t)  *  h(t  +  t)  dt .  Eq.  3.1.1 

The  correlation  value  depends  on  the  delay  value  i.  For  G(f),  the  Fourier 
transform  of  g(t)  and  H(f),  the  Fourier  transform  of  h(t),  then 

Corr(g,h)  =  F"1  (G(f)  H*  (f))  Eq.  3.1.2 

where  F-1  is  the  inverse  Fourier  transform  and  the  asterisk  refers  to  the 
complex  conjugate  of  the  transform  value.  In  the  implementation,  the 
mother  wavelets  were  pre-computed  for  various  alpha  (scale)  values  at 
zero  delay.  The  resulting  functions  were  transformed  and  stored  as 
transformed  values.  The  inner  product  values  at  each  scale  and  delay 
were  performed  by  computing  the  transform  of  the  input  function  (g)  once 
and  then  performing  the  complex  product  operations  for  each  alpha 
representation.  The  inverse  transforms  of  the  products  produce  the 
transform  values  as  a  function  of  both  scale  and  delay. 

As  an  example  of  the  processing  change,  a  512  by  512  pixel  image  processed  using 
the  WAV_MAC  algorithm  on  a  Silicon  Graghics  IP12  Processor  (36MHz), 
required  59  hours  processing  time  before  implementation  of  the  FFT  method  and 
5.4  hours  after  implementation  of  the  FFT  method. 

3.2  Orthonormal  Wavelets 

In  addition  to  the  general  frame  wavelets,  the  orthonormal  4  point  Daubechies 
wavelet  was  also  implemented  and  is  used  exclusively  in  the  WAV_TEN, 
WAV_SOF,  and  WAV_VOL  tools  and  implemented  with  general  frame  wavelets 
in  the  WAV_CEP  tool.  The  pyramidal  form  of  the  transform  and  its  inverse 
produce  computation  speeds  that  approach  that  of  FFTs  and  permits  fast 
inspection  of  the  effects  of  the  tool  being  used  on  the  input  image.  Until  such  time 
that  the  general  frame  method  can  be  sped  up  to  yield  comparable  processing 
times,  the  orthonormal  wavelets  are  the  method  of  choice. 
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4.0  WAVELET  TOOLS 


4.1  WAV_MAC 

The  purpose  of  this  tool  is  to  improve  visualization  by  imparting  a 
pseudo-coloration  of  the  image.  The  colors  are  derived  from  the  absolute  value  of 
the  wavelet  transform  of  the  image  partitioned  into  groups  formed  by  consecutive 
scales.  The  lower  consecutive  scales  are  assigned  red,  the  mid  range  consecutive 
scales  are  assigned  green,  and  the  high  consecutive  scales  are  assigned  blue. 

The  color  weights  at  each  pixel  sum  to  one  and  are  used  as  multipliers  of  the 
original  image  intensity  values.  Thus,  features  in  the  image  are  colorized  by 
scale.  With  clipping,  shadow  areas  remain  black  and  bright  targets  are 
rendered  as  white. 

WAVJMAC  utilizes  the  FFT  method  (described  in  the  Wavelets  Functions  and 
Implementation  section)  to  decrease  processing  time.  With  this  method  there  are 
8  (MAX_ALPHA)  scales  produced  in  the  wavelet  domain.  Scales  0-2  are 
assigned  to  RED,  scales  3-5  are  assigned  GREEN,  and  scales  6-7  are  assigned 
BLUE.  These  scale  boundaries  are  set  by  the  "#define"  statements  in  the 
IPWavMACO  function  source  code  and  can  be  changed  by  the  user  for  varying 
results.  A  block  diagram  and  example  of  WAVJVIAC  are  given  in  Figure  4.1-1. 


4.2  WAVJTEN 

The  purpose  of  this  tool  is  to  demonstrate  the  potential  for  image  compression 
using  wavelets.  The  tool  employs  a  two  dimensional  four  point  Daubechies 
wavelet  routine  and  a  histograming  technique  based  on  the  log  of  the  absolute 
value  of  the  wavelet  transform  data  to  choose  a  threshold  level.  A  percentage  of 
the  data  to  be  retained  is  input  to  the  routine  and  this  value  is  translated  to  a 
threshold  value  (transform  data  below  this  level  is  set  to  zero).  The  effects  on  the 
image  are  immediately  visualized  by  reconstruction  of  the  image  with  the 
thresholded  transform  data.  The  transform  is  output  for  storage  to  an  external 
file  so  that  subsequent  processing  can  be  done.  The  data  are  stored  as  floating 
point  values.  It  is  noted  that  even  though  the  input  data  is  8  bits  (size  of  a 
character,  1  byte),  a  larger  dynamic  range  is  required  to  represent  the 
transformed  data  (size  of  a  float,  4bytes).  Studies  would  have  to  be  conducted  to 
determine  the  word  size  to  represent  the  data.  Further  enhancements  of  this  tool 
would  be  the  incorporation  of  the  compression  algorithms  operation  on  the  sparse 
data  set.  A  block  diagram  and  example  of  WAVTEN  are  given  in  Figure  4.2-1. 
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Figure  4.2-1  WAV_TEN  (Wavelet  Transform  Encoding) 


4.3  WAV_SOF 


The  Mallat  decomposition  provides  a  means  for  visualization  of  local  horizontal, 
vertical  or  a  combination  of  these  data  in  an  image.  The  process  utilizes  a  four 
point  Daubechies  wavelet  transform  operating  on  the  image  and  a  transposed 
version  of  the  image.  The  transform  representation  provides  half  of  the  total 
number  of  points  (assumed  a  power  of  2)  to  represent  the  smallest  scale,  half 
again  for  the  next  larger  scale  and  so  on.  The  image  is  processed  with  zero  fill 
elements  such  that  the  picture  fills  a  matrix  to  the  nearest  power  of  2  in  each 
dimension.  For  a  matrix  size  2m  by  2n,  the  transform  is  presented  in  four 
quadrants  of  size  2(m_1)  by  2(n_1).  By  operation  on  the  image  and  its  transpose, 
one  quadrant  represents  horizontal  data,  another  vertical  data,  a  third  is  a 
combination  of  horizontal  and  vertical,  and  the  fourth  a  replication  of  the  original 
image  but  reduced  in  resolution  in  each  dimension  by  two.  A  replication  occurs 
because  the  transform  is  orthonormal  and  applying  the  transform  to  the 
transpose  of  the  transform  results  in  an  identity.  The  presentation  shows  the 
absolute  value  of  the  transform  which  is  scaled  to  fill  the  0-255  grey-scale  range. 
Zero  fill  regions  in  the  presentation  (represented  by  black  bands)  could  be 
eliminated  in  future  versions  of  the  IPWavSOFO  function. 

The  WAV_SOF  program  is  designed  to  perform  repeated  Daubechies  wavelet 
transforms  on  the  reduced  resolution  replicate  of  the  original  image  present  in 
the  transform  output.  The  number  of  times  the  transform  is  repeated  is 
controlled  by  argument  "k_max”  of  the  function  IPWavSOFO.  Testing  of  the 
function  indicate  that  "k_max"  >  3  gives  limited  to  no  added  utility.  The 
argument  "k_max"  is  presently  hardcoded  to  1  in  the  IPToolKit  file 
"enhance_handler.c".  It  is  suggested  that  if  future  changes  are  made  to  the 
IPToolKit  menus  to  allow  multiple  selections  of  "k_max",  the  values  be  limited  to 
1  through  3.  A  block  diagram  and  example  of  WAV_SOF  are  given  in  Figure 
4.3-1. 


4.4  WAVjCEP 

The  wavelet  transforms  offer  a  method  of  enhancing  some  image  features  on  the 
basis  of  object  size  (or  scale).  Histograms  of  transform  energy  (the  square  of  the 
transform  values  accumulated  over  delay  for  a  particular  scale)  versus  scale 
value  exhibit  a  general  fall  off  from  large  to  small  scale  values  not  unlike  Power 
Spectral  Densities  (PSD)  derived  in  the  Fourier  domain.  Logarithmic 
compression  of  the  transform  values  (using  the  Cepstrum  calculation)  reduces 
the  relative  ratio  of  the  large  to  small  scale  data.  The  resultant  effect  is 
enhancement  of  small  features  and  edges. 
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Figure  4.3-1  WAV_SOF  (Wavelet  Scale-Orientation  Filter) 


The  dyadic  grid  technique  described  earlier  was  implemented  for  the  general 
frame  wavelets  of  IPWavCepstrum().  Though  this  technique  decreased  execution 
time,  the  Morlet  wavelet  waveform  could  not  be  accurately  implemented.  Thus, 
the  results  for  the  Morlet  wavelet  show  little  utility.  A  block  diagram  and  example 
of  WAV_CEP  are  given  in  Figure  4.4-1. 


4.5  WAV_SCA 

The  wavelet  transforms  are  labeled  according  to  delay  (the  location  of  a  local  inner 
product  calculation)  and  scale  (the  extent  of  the  wavelet  function  for  a  compact 
wavelet  and  the  region  containing  most  of  the  wavelet  function  energy  for  the 
general  form  case). 

A  single  image  vector  derived  from  a  single  row  in  the  image  is  processed  in  this 
routine  over  a  delay  range  corresponding  to  the  vector  length  and  over  a  range  of 
scales.  A  plot  of  transform  output  at  a  particular  scale  versus  delay  produces  a 
scalogram.  Certain  signal  aspects  such  as  bright  point  objects  exhibit  local  high 
amplitude  concentrations  in  the  scalogram. 

Because  of  the  large'  absolute  amplitude  dynamic  range  that  occurs  with  scale 
change,  a  non-linear  mapping  and  color  are  required  to  maximize  the 
visualization  of  the  data.  Incorporation  of  the  required  algorithms  are  reserved 
for  a  future  update.  The  current  process  produces  a  quantized  gray  scale 
presentation  using  function  "sca_pix_hist".  A  block  diagram  and  example  of 
WAV_SCA  are  given  in  Figure  4.5-1. 


4.6  WAV  VOL 

The  wavelet  transform  is  parameterized  (labeled)  according  to  amplitude,  scale 
and  delay.  For  some  simple  signal  structures  such  as  a  range  profile  output  of  a 
High  Range  Resolution  radar,  the  transform  over  some  delay  range  (representing 
the  profile  length  relative  to  a  specific  starting  point),  over  some  scales,  and  above 
some  threshold  level  determine  a  volume  in  transform  space.  With  a  judicious 
choice  of  parameters,  the  inverse  transform  of  the  selected  transform  values 
exhibit  a  signature  with  less  noise  than  the  original  and,  possibly,  somewhat  less 
resolution.  For  more  complex  structures  such  as  a  2D  image,  the  delay  is  not  as 
obvious  and  may  be  represented  by  rectangular  regions  (windows)  to 
accommodate  many  objects  of  similar  size  dispersed  over  the  image.  In  the 
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Scalogram 


current  volume  pass  implementation,  the  selection  of  particular  windows  is 
reserved  as  a  post-processing  function  (and  is  currently  not  implemented).  The 
central  focus  of  the  processing  is  to  manipulate  the  transform  according  to  the 
scale  data.  Controls  are  provided  to  select  scales  in  forms  similar  to  that  in  the 
frequency  domain:  low  pass  (largest  to  moderate  scales  in  consecutive  order), 
high  pass  (moderate  to  smallest  scale  in  consecutive  order),  band  pass  (a  selected 
set  with  zero  values  for  scale  data  not  in  the  set)  and  stop  band  (a  selected  set  of 
scales  set  to  zero  with  others  left  intact).  An  additional  control  is  provided  so  that 
the  above  selections  can  be  added  or  subtracted  from  the  total  transform 
corresponding,  for  example,  to  boosting  the  data  amplitude  over  some  region  in 
the  frequency  spectrum.  These  operations  produce  images  with  certain  objects 
and  features  which  are  accentuated  on  the  basis  of  scale  (or  object  size)  and 
amplitude.  Further  operations  on  the  enhanced  image  could  by  provided  by 
window  filters  as  indicated  above.  A  block  diagram  and  example  of  WAV_VOL 
are  given  in  Figure  4.6-1. 
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Figure  4.6-1  WAV_VOL  (Wavelet  Volumepass  Filtering) 


5.0  SUMMARY 


This  report  has  presented  background  information  and  some  of  the  tool  specifics 
on  the  application  of  wavelet  transform  methods  in  the  enhancement  and 
visualization  of  digital  reconnaissance  imagery.  The  tools  exploit  the  localization 
properties  of  the  transform  and  the  ability  to  manipulate  data  based  on  scale  (or 
object  size).  The  WAV_TEN  tool,  for  example,  exploits  the  sparse  nature  of  the 
discrete  wavelet  transform  and  the  ability  to  represent  the  transform  with  a  small 
subset  of  the  computed  values  (the  other  values  being  set  to  zero).  Other  tools  show 
the  effect  of  the  mother  wavelet  on  the  processed  image  by  way  of  the  presentation 
of  specific  image  features. 


6.0  RECOMMENDATIONS 

Methods  should  be  sought  to  improve  the  throughput  rates  for  the  wavelets  which 
are  currently  being  computed  using  general  frame  methods  (and  some  adjuncts 
to  these  methods).  Further  study  is  required  to  determine  the  construction  and 
approximation  techniques  for  these  waveforms  that  may  be  applicable. 

Some  minor  cosmetic  improvements  in  presentation  of  the  data  could  be  made  on 
the  IPWavSOFO  function  to  account  for  image  sizes  which  are  not  powers  of  two 
(the  black  regions  resulting  from  zero  fill  of  the  original  data  can  be  eliminated  by 
computation  of  their  location  in  transform  space).  Similarly,  non-linear  intensity 
mapping  and  color  would  enhance  the  IPWavScalogramO  function  presentation. 

The  WAV_MAC  tool  is  hardcoded  to  divide  the  wavelet  domain  scales  into  3 
approximately  equal  groups  (for  red,  green  and  blue).  Future  upgrade  could 
include  an  interactive  scroll  bar  to  vary  the  scales  into  unequal  groups. 
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